from math import *
import matplotlib.pyplot as plt
from matplotlib import font_manager as fm
from matplotlib import rcParams
from matplotlib.ticker import MultipleLocator


def XY_Axis(x_start, x_end, y_start, y_end):
    plt.rcParams['figure.dpi'] = 1200

    fig = plt.figure(figsize=(5, 4))
    ax = fig.add_subplot(111)

    xmajor = MultipleLocator(0.5 * x_end)
    ymajor_1 = MultipleLocator(y_end * 0.25)

    ax.xaxis.set_major_locator(xmajor)
    ax.tick_params(axis="x", direction="out", which="major", labelsize=9, length=3)
    plt.xticks([x_start, (x_start + x_end)/2, x_end], ['0', 'arcsin($r/L$)', '2arcsin($r/L$)'])

    ax.set_yticks([])

    plt.text(x_end * 0.488, y_start - 0.36, "θ", fontproperties = fm.FontProperties(fname = 'C:/Windows/Fonts/TIMESI.TTF', size = 14, style = 'italic'))
    plt.text(x_start - 0.2, y_end * 0.465, "f$_{4}$(θ)", fontproperties = fm.FontProperties(fname = 'C:/Windows/Fonts/TIMESI.TTF', size = 14, style = 'italic'), rotation = 90)


def ComputingDenominator(L, r, mAngle, N1, N2):
    L2 = L * L
    r2 = r * r
    r4 = r2 * r2
    r5 = r4 * r
    r6 = r4 * r2
    A = mAngle / N1
    iDenominator = 0.0
    for i in range(0, N1):
        a = (i + 0.5) * A
        s = sin(a)
        s2 = s * s
        c = cos(a)
        c2 = c * c
        sc = s * c
        t = tan(a)
        Ls = L * s
        Lc = L * c
        Lsc = Ls * c
        hLc2 = 0.5 * L * c2
        hLs = 0.5 * Ls
        hLsc = 0.5 * Lsc
        mLength = -0.5 * L * abs(t) + r / c
        B = mLength / N2
        jDenominator = 0.0
        for j in range(0, N2):
            b = (j + 0.5) * B
            bc = b * c
            bc2 = b * c2
            bsc = b * sc
            Et1 = hLs - bc
            Delta1 = r2 - Et1 * Et1
            if Delta1 < 0.0:
                print("i = {}, j = {}, Delta1 = {}".format(i, j, Delta1))
                Esqrt1 = 0.0
            else:
                Esqrt1 = sqrt(Delta1)
            x1 = -bsc - hLc2 - c * Esqrt1
            y1 = bc2 - hLsc - s * Esqrt1
            Et2 = hLs + bc
            Delta2 = r2 - Et2 * Et2
            if Delta2 < 0.0:
                print("i = {}, j = {}, Delta2 = {}".format(i, j, Delta2))
                Esqrt2 = 0.0
            else:
                Esqrt2 = sqrt(Delta2)
            x3 = -bsc + hLc2 - c * Esqrt2
            y3 = bc2 + hLsc - s * Esqrt2

            ta1 = Ls - 2 * bc
            ta2 = ta1 * ta1
            Delta3 = 4 * r2 - ta2
            if Delta3 < 0.0:
                print("i = {}, j = {}, Delta3 = {}".format(i, j, Delta3))
                t1 = 0.0
            else:
                t1 = sqrt(Delta3)
            x1c = x1 * c
            y1s = y1 * s
            x1cPy1s = x1c + y1s
            x12Py12Pr2 = x1 * x1 + y1 * y1 + r2
            r4E1 = 0.2
            r4E2 = x1cPy1s
            r4E3 = (4 * x1cPy1s * x1cPy1s + 2 * x12Py12Pr2) / 3
            r4E4 = 2 * x1cPy1s * x12Py12Pr2
            r4E5 = x12Py12Pr2 * x12Py12Pr2
            
            tb1 = Ls + 2 * bc
            tb2 = tb1 * tb1
            Delta4 = 4 * r2 - tb2
            if Delta4 < 0.0:
                print("i = {}, j = {}, Delta4 = {}".format(i, j, Delta4))
                t2 = 0.0
            else:
                t2 = sqrt(Delta4)
            x3c = x3 * c
            y3s = y3 * s
            x3cPy3s = x3c + y3s
            x32Py32Pr2 = x3 * x3 + y3 * y3 + r2
            r4E6 = x3cPy3s
            r4E7 = (4 * x3cPy3s * x3cPy3s + 2 * x32Py32Pr2) / 3
            r4E8 = 2 * x3cPy3s * x32Py32Pr2
            r4E9 = x32Py32Pr2 * x32Py32Pr2

            r4w1d1 = r4E1 * (t1 ** 5) + r4E2 * (t1 ** 4) + r4E3 * (t1 ** 3) \
                     + r4E4 * (t1 ** 2) + r4E5 * t1
            r4w2d2 = r4E1 * (t2 ** 5) + r4E6 * (t2 ** 4) + r4E7 * (t2 ** 3) \
                     + r4E8 * (t2 ** 2) + r4E9 * t2
            r8w12d12 = r4w1d1 * r4w2d2
            jDenominator += r8w12d12
        jDenominator = jDenominator * B / (r ** 8)
        iDenominator += jDenominator
    Denominator = iDenominator * A * 2.0
    return(Denominator)


def D_PDF_TwoP(L, r, mAngle, N1, N2):
    Denominator = ComputingDenominator(L, r, mAngle, N1, N2)
    L2 = L * L
    r2 = r * r
    r4 = r2 * r2
    r5 = r4 * r
    r6 = r4 * r2
    A = mAngle / N1
    iNumerator = 0.0
    angleArray = []
    Function = []
    Fmax = 0.0
    for i in range(-N1, N1):
        a = (i + 0.5) * A
        angleArray.append(a + mAngle)
        s = sin(a)
        s2 = s * s
        c = cos(a)
        c2 = c * c
        sc = s * c
        t = tan(a)
        Ls = L * s
        Lc = L * c
        Lsc = Ls * c
        hLc2 = 0.5 * L * c2
        hLs = 0.5 * Ls
        hLsc = 0.5 * Lsc
        mLength = -0.5 * L * abs(t) + r / c
        B = mLength / N2
        jNumerator = 0.0
        for j in range(0, N2):
            b = (j + 0.5) * B
            bc = b * c
            bc2 = b * c2
            bsc = b * sc
            Et1 = hLs - bc
            Esqrt1 = sqrt(r2 - Et1 * Et1)
            x1 = -bsc - hLc2 - c * Esqrt1
            y1 = bc2 - hLsc - s * Esqrt1
            Et2 = hLs + bc
            Esqrt2 = sqrt(r2 - Et2 * Et2)
            x3 = -bsc + hLc2 - c * Esqrt2
            y3 = bc2 + hLsc - s * Esqrt2

            ta1 = Ls - 2 * bc
            ta2 = ta1 * ta1
            t1 = sqrt(4 * r2 - ta2)
            x1c = x1 * c
            y1s = y1 * s
            x1cPy1s = x1c + y1s
            x12Py12Pr2 = x1 * x1 + y1 * y1 + r2
            r4E1 = 0.2
            r4E2 = x1cPy1s
            r4E3 = (4 * x1cPy1s * x1cPy1s + 2 * x12Py12Pr2) / 3
            r4E4 = 2 * x1cPy1s * x12Py12Pr2
            r4E5 = x12Py12Pr2 * x12Py12Pr2
            
            tb1 = Ls + 2 * bc
            tb2 = tb1 * tb1
            t2 = sqrt(4 * r2 - tb2)
            x3c = x3 * c
            y3s = y3 * s
            x3cPy3s = x3c + y3s
            x32Py32Pr2 = x3 * x3 + y3 * y3 + r2
            r4E6 = x3cPy3s
            r4E7 = (4 * x3cPy3s * x3cPy3s + 2 * x32Py32Pr2) / 3
            r4E8 = 2 * x3cPy3s * x32Py32Pr2
            r4E9 = x32Py32Pr2 * x32Py32Pr2

            r4w1d1 = r4E1 * (t1 ** 5) + r4E2 * (t1 ** 4) + r4E3 * (t1 ** 3) \
                     + r4E4 * (t1 ** 2) + r4E5 * t1
            r4w2d2 = r4E1 * (t2 ** 5) + r4E6 * (t2 ** 4) + r4E7 * (t2 ** 3) \
                     + r4E8 * (t2 ** 2) + r4E9 * t2
            r8w12d12 = r4w1d1 * r4w2d2
            jNumerator += r8w12d12
        iNumerator = jNumerator * B / (r ** 8)
        tFunction = iNumerator / Denominator
        Function.append(tFunction)
        if Fmax < tFunction:
            Fmax = tFunction
    XY_Axis(0.0, 2 * mAngle, 0.0, Fmax)
    plt.plot(angleArray, Function, "k-", linewidth=1, linestyle="-")
    plt.savefig("Fig6.jpg")
    plt.show()


L = 100
r = 35
mAngle = asin(2 * r / L)
N1 = 1000
N2 = 1000
D_PDF_TwoP(L, r, mAngle, N1, N2)
